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The propagator for the activity in a broad class of self-organized critical models 
obeys an imaginary-time Schrodinger equation with a nonlocal, history-dependent 
potential representing memory. Consequently, the probability for an avalanche to 
spread beyond a distance r in time t has an anomalous tail exp 

\_ C X 1/{D-1)] for 

x = r D /t ^> 1 and D > 2, indicative of glassy dynamics. The theory is verified for an 
exactly solvable model, where D = 4 and C = 3/4, and for the Bak-Sneppen model 
where it is tested numerically. 
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At the simplest level, the dynamics of large, extended systems may be either fast and 
chaotic, or slow and possibly frozen. In the first case information about the current micro- 
scopic state decreases rapidly in time; in the second the system acquires long-term memory. 
As some tuning parameter is varied, there could generally be a critical transition between 
these two phases where the dynamics takes place intermittently in terms of avalanches. Pos- 
sible examples include depinning transitions Jy], dynamic glass transitions 0, and certain 
transitions in cellular automata |§. Many slowly driven, dissipative systems may naturally 
operate at this transition; the criticality is self-organized (SOC) j§]. Such systems evolve in 
terms of scale-free avalanches sweeping the system. This picture can be verified in a broad 
class of SOC models [|lj including interface || and flux-line depinning ||, and certain aspects 
of biological evolution [7||§. These models represent constant- velocity freezing transitions 
where the velocity v — > of motion is the order parameter. The intermittent, punctuated 
equilibrium behavior at the phase transition (v = 0) combines features of frozen systems 
with those of chaotic ones. Such systems can remember the past because of the long periods 
of stasis allowing them to preserve what has been learned through history; they can evolve 
rapidly and adapt because of the intermittent bursts of activity, which span all length scales 
up to the system size. 

The freezing transition is an attractor for the dynamics when the site which is updated 
at each instant has the current, globally extremal value of some force. The extremal update 
causes that site as well as neighboring sites to get new forces according to specific model- 
dependent rules. In the stationary critical state, the pattern of activity is correlated and 
can be described in terms of scale free avalanches. The avalanche propagator G(r, t) is 
the conditional probability to find the updating activity at r at time t given a surviving 
avalanche which started at the origin at time t = 0. Knowledge of the propagator as well as 
the distribution of avalanche durations contains all of the information on the macroscopic 
behavior of the system. Any universal quantities which may exists would appear in the 
asymptotic long length and long time limits of these functions. 

Our main result is that the asymptotic behavior of the propagator, G, is governed by an 
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imaginary-time Schrodinger equation with a nonlocal, history-dependent potential 

^M = Vr 2 G(r,t)+ jW(t-t')G(r,t')dt' . (1) 

As the system builds up to the critical state the potential becomes scale free, 

V(t) ~ r a , (t -> oo). (2) 

Such dynamics always leads to subdiffusive spreading of the avalanche, with an avalanche 
dimension D > 2, where D — 2/ (a — 1). The equation governing the asymptotic behavior 

of the propagator adheres to an anomalous non-Gaussian diffusion form, 
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G(r,t) ~ A(r,t)exp 
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where C is a constant and A(r, t) contains only powers of r and t. For the Bak-Sneppen 
evolution model we have verified this prediction numerically as shown in Figs. 1 and 2. 
In a similar model, the multi-trait evolution model [Q, Eqs. (|1]-§|) can be explicitly derived 
from microscopic rules. In this case, the nonlocal potential V(t) is obtained exactly; it is 
the distribution of avalanche durations. The resulting equation for the propagator is solved 
exactly with a = 3/2 to give Eq. (g) with D = A and C = 3/4. 

Similar dynamical equations with memory describing aging of glasses and spin glasses 
within mode-coupling theory have been introduced recently fHU]. This similarity suggests a 
possible connection between ergodicity breaking in glasses and the punctuated equilibrium 
behavior of self-organized critical systems. Perhaps, like SOC, the dynamics of glasses is at 
the edge separating frozen and chaotic phases, although previous results testing this view 
have been mixed [jK],|ll|] . 



The simplest extremal model is the Bak-Sneppen evolution model defined as follows [0 : 
Put a random number at each site on a d dimensional lattice. At each time step, replace 
the smallest random number in the system and its 2d neighbors with new random numbers 
drawn independently from a flat distribution between and 1. After a long transient the 
model self-organizes to a critical state where almost all sites have values above a threshold 



A c . An avalanche of duration t begins at time when the minimal random number in the 
system equals A c and ends at the first instant in time, t, when the minimal random number 
reaches A c again. The sequence of values below A c form bursts which span all spatial and 
temporal extents. Reducing the threshold level A < A c divides the critical avalanches into 
smaller A subavalanches (with a cutoff) which form a hierarchy For any A < Ao < A c 
the A subavalanche must persist longer than the A sub-avalanche that it contains. In all 
extremal models, the avalanche hierarchy is characterized by certain power laws Thus, 
a mechanism for long-term memory is apparent. 

We propose that the propagator G for the avalanches in these models obeys Eqs. 
for r D ^> t ^> 1. Current activity has some probability to activate neighboring sites at 
subsequent nearby times. In an isotropic system, this by itself would give ordinary diffusion 
as represented by Eq. ([!]) with V = 0. The memory term with nonzero V explicitly 
breaks invariance with respect to shifts in time. The history-dependent potential V(t) is 
model dependent. This potential accounts for the avalanche hierarchy described above. 
We conjecture that in all extremal SOC models V(t) is a power law for large arguments 
t; otherwise there would be a characteristic scale where the dynamics would no longer 
be dominated by memory, and where the avalanche hierarchy would break down. The 
propagator at time t gets contributions from all previous times (t f < t) weighted by V(t — t'), 
which is scale free. Since we consider only the behavior in the tail of the distribution, 
r D /t ^> 1, an equation linear in G is sufficient. A complete equation for the propagator, 
such that probability is conserved, would be far more complicated and would most likely be 
nonlinear [as in Eq. (|^) for the explicit example below]. 

We checked Eq. (|]) numerically for the two dimensional Bak-Sneppen model using A c = 
0.328855 and D = 2.92 as determined in Ref. M. Based on previous numerical studies of 



the multi-trait model [I2"|, we assume that G(r,t) has the scaling form |13 



G(r, 0-^(^1 (r D /t»l). (4) 



v t t 

At each update t we determine the location r of the minimal random number in a surviving 



avalanche that started at t = 0, r = 0. We bin counts as a function of log 2 1 and log 2 (r D /t) 
for avalanches of duration t < 10 9 ; see Fig. [I]. The data indicates that the counting rates 
for each bin rise with t s , 6 ~ 0.7, for each value of r D /t. We then "collapse" the data to 
determine g(x) ~ £ 1-a (j(r, t) with the scaling variable x = r D /t by averaging the data for 
each x. We attempted to numerically determine the asymptotic tail of the function g for 
large values of x by assuming in accordance with Eq. (||) that 



g{x) ~ exp 
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x>l. (5) 



Since the estimate of D' depends on taking two logarithms of g to extract the asymptotic 
behavior for large x, it is not very accurate. From the extrapolation in Fig. we estimate 
that 1/(D' — 1) = 0.55 ± 0.10, consistent with the expected result of 0.52 obtained for 
D = 2.92 [|IJ. Note that the apparent convergence of the extrapolation justifies Eq. fl5|). 

Eqs. (|iy^) can be derived by considering a similar model which is analytically tractable, 
the multi-trait evolution model. This model has all of the characteristic punctuated equi- 
librium behavior of SOC including dimension- dependent critical exponents. Now, each site 
on the lattice has a collection of M random numbers. The smallest among all M x L d is 
chosen. That particular random number and one random number out of M at each of the 
2d neighboring sites are replaced by new random numbers independently drawn from a flat 
distribution in the unit interval. When M — 1, the Bak-Sneppen model is recovered. When 
M — > oo the equation of motion can be obtained explicitly and solved. In what follows we 
choose M = oo, d = 1, and modify the model slightly by replacing the minimal random 
number with the value 1. Then, A c = 1/2 exactly ||. 

Let F(r,t) be the probability for a A c avalanche to survive precisely t steps and to have 
affected a particular site of distance r from its origin. Conceptually, the quantity F(r, t) may 
roughly correspond to an envelope function of the propagator G(r,t). Due to the lack of 
any scale in the model, it is plausible that the asymptotic behavior of G and F is identical, 
as comparison with numerical calculations suggests |T2| . 



To simplify the derivation of an equation for F, we consider P(r,t), the probability that 



the A c avalanche dies precisely after t updates and does not affect a particular site r away 
from the origin of the avalanche. Clearly, 

P(r,t) = P(r = oo,t) - F(r,t), (6) 

where P(r = oo,t) = P(t). Since no avalanche of finite duration spreads to an infinitely 
distant site, P(t) is the probability for an avalanche to end exactly after t updates, or its 
lifetime distribution. In Ref. || it was shown that in this model P(t) is given by the return 
time distribution for a random walk, 

p w - vlmU ~^ 1/2r3/2 for f>>1 • (7 » 

Since we consider the avalanche to start with a single active value at r = and t = 0, 
P{r = 0,t) =0 for all t > 0, and P(r,t = 0) = for all r. The remaining properties of 
a A c avalanche can be deduced from the properties of avalanches that ensue after the first 
update. It will terminate with probability (1 — A c ) 2 = 1/4 after the first update. Thus, 
P(r, t = 1) = 1/4. For avalanches surviving until t > 2 we find for r > 1 

P(r, t) = - [P(r - 1, t - 1) + P(r + 1, t - 1)] 

+^F(r-M')P(r + l,f-l-t') (8) 

4 t'=0 

expressing that the first update may create exactly one new active barrier with probability 
A c (l — A c ) = 1/4 either to the left or to the right of the origin, and also that the first update 
may create two new active barriers with probability A 2 = 1/4 to the left and the right of 
the origin. Then, the properties of the original avalanche of duration t are related to the 
properties of all combinations of two avalanches of combined duration t — 1. Both of these 
avalanches evolve in a statistically independent manner for M = oo. For any such pair, 
the probability to not affect the chosen site of distance r from the origin is given simply 
by the product of the probabilities for the two ensuing avalanches to not affect a chosen 
site of distance r — 1 or r + 1, respectively. In this respect, the theory is similar to the 
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mode-coupling theory used in Refs. 0||, except that in our case the theory is exact. Note 
that setting r = oo gives a nonlinear equation for P(t) whose solution is Eq. ([T|). 

The equation governing the envelope function F(r, t) is obtained by inserting Eq. (||) 
into Eq. (§). This gives F{0,t) = P(t), F(r,0) = 0, F(r > l,t = 1) = 0, and for t > 1, 
r > 1, 

F(r, t + 1) = i [F(r -l,t) + F(r + 1, t)] 

+]i: p (t-t') [F(r-l,t') + F(r+l,t')] 

4 t'=0 

-^F(r-l,^(r + M-t'). (9) 

4 t'=0 

One can show that F(r,t) — ► for r D ^> t ^> 1, D = 4, sufficiently fast such that the 
nonlinear term in Eq. (§) can be neglected. We can take the continuum limit and obtain 
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dF{r, t) 



dt ^ 2 rF(r, t) + \ jT V(t - t')F(r, t')dt', (10) 

with a nonlocal kernel V(t) = P(t) — 28(t). Aside from trivial constants, Eq. (0) is identical 
to Eq. (D). In the absence of the history-dependent potential, the system would be purely 
diffusive with a Gaussian tail F ~ e _r-2//2 *. In its presence, the probability to have reached a 
site at distance r at time t gets contributions from avalanches that reached r at earlier times 
if < t weighted according to P(t — if) which has a power-law tail. Avalanches contributing 
to F(r, t) consist of sub- avalanches, one of which reaches r in time t' while the combined 
duration of the other is t — t'. The avalanche hierarchy gives a hierarchy of time scales 
changing the relaxation dynamics to be non-Gaussian. 

Using a Laplace transform, G(r,y) = / °° dt e~ yt G(r,t), Eq. (|I|) turns into an ordinary 
second-order differential equation in r, 

V?<2(r,y)~ [2y-V(y)]G(r,y), (11) 

where V(y) is the Laplace transform of V(t); V(y) ~ y^ 1 for small y. In our particular 
model where V(t) ~ P(t) ~ t -3 / 2 , a = 3/2. We require that a > 1 so that V can be 



normalized. On the other hand, we require a < 2, so that V(y) ^> y for y — > (i.e. for 
t — > oo) in Eq. (pTT|) . Then in the limit of large times, the effect of the history-dependent 
potential dominates over the time derivative in Eq. ([!]), signaling the deviation from simple 
diffusive behavior. 

The solution of Eq. (11) which decreases for large r is 

~ r a—l~ 

G(r,y) ~ exp — Ary^~ , (12) 

where A is a constant, and we assume that any prefactor in Eq. (0) is a sufficiently well- 
behaved function near y = 0. 

The inverse Laplace transform yields a contour integral with a contour extending just to 
the right (77 > 0) of the imaginary y-axis: 

G(r,t) ~ / —exp yt-Ary— . (13) 



For r/t^ 1 ^ 2 — > 00 the integral in Eq. (|T3|) has a nontrivial saddle point such that G(r,t) 
is exponentially cut off, determining the avalanche dimension D = —^-r. With 1 < a < 2 
we have D > 2, i. e. the avalanche dynamics as described by Eq. ([I]) always evolves in a 
subdiffusive manner. We note that almost all isotropic, extremal SOC models have been 
found to have D > 2 ||. The notable exception is the one-dimensional interface depinning 
model proposed by Sneppen [Q, which therefore cannot be described by this theory. In our 
particular case we find D = 4 from a = 3/2. In the limit r D /t — > 00, we can perform 
a steepest-descent analysis of the integral Jl5| and obtain the non-Gaussian tail given in 
Eq. (|). In the same way, the complete leading asymptotic behavior for F in Eq. (9) can 
be calculated explicitly to give [|i~2] 



F M) ~ VV rl (j)' 6 (r 4 >t>l) • (14) 

Note that previous continuum theories for SOC, in particular the theory of singular 
diffusion by Carlson and others, specifically apply to systems with a conserved quantity; 



whereas in the cases considered here, there is no conservation law. This does not rule out 
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the possibility that equations similar to Eq. (1-3) could describe avalanches in, for instance, 
sand pile models, which involve transport of a conserved quantity. Critical exponents for 



avalanches in a limited local sand pile model have been related to the singularity in the 
diffusion equation |TJ| where the avalanches are averaged in a coarse graining procedure. 
This topic is currently under investigation. 

It has previously been speculated that glassy dynamics may take place near a (self- 
organized) critical point ||10|| , with mixed results flUf . Here we find an anomalous tail in the 
probability distribution for the activity in SOC which is characteristic of glassy systems. 
For example, the directed polymer in a random media, which contains many features of 
frustrated random systems, also has a non-Gaussian tail for G(x,t) [pl|. Our problem, 
however, is inherently dynamical with no quenched disorder. From a different viewpoint, 
Stein and Newman |2(J have put forward a picture of slow dynamics on a high dimensional 
rugged fitness landscape based on an invasion percolation or SOC picture. Finally, the 
mode-coupling equations for glasses H such as those derived for the non-ergodic dynamics 
of the random phase Sine-Gordon model || are similar to what we find for avalanches in 
self-organized criticality. It would be interesting to see if the mode coupling equations in 
glasses could be interpreted in terms of avalanches sweeping through a critical system. 
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FIGURES 

FIG. 1. Counting rates for the activity plotted as a function of time t for the two dimensional 
Bak-Sneppen model. Each curve is labeled in descending order by i = |~log 2 x\ where x = r D /t is 
the scaling variable [see Eq. @]. 
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FIG. 2. Extrapolation of the data in Fig. |T| according to Eq. (|5|). We plot for each i = |~log 2 x] 
the projected value of 1/(D' — 1) ~ In [— lng(x)] / In x as a function of 1/lnx and determine the 
value for 1/(D' — 1) = 0.55 ± 0.10 from the extrapolation of this sequence for x — > oo, i. e. at 
1/lnx = 0. Error bars originate from the averaging of the data in Fig. |l| for each i. 
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